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Abstract  A Large  Eddy  Simulation  (LES)  methodology  has  been  developed  for 
supersonic  turbulent  flows  with  strong  shock  boundary  layer  interaction. 
Results  are  presented  for  an  expansion-compression  corner  at  Mach  3 
and  compared  with  experimental  data. 


Introduction 

The  interaction  of  shock  waves  and  turbulent  boundary  layers  is  a 
common  and  important  phenomenon  in  aerodynamics,  and  has  been 
studied  extensively  (Settles  and  Dolling,  1990;  Zheltovodov,  1996).  Con- 
ventional Reynolds- averaged  Navier-Stokes  methods  have  been  unable  to 
accurately  predict  separated  shock  wave  turbulent  boundary  layer  inter- 
actions (Knight  and  Degrez,  1998).  Recently,  Large  Eddy  Simulation 
(LES)  and  Direct  Numerical  Simulation  (DNS)  have  been  applied  to 
shock  wave  turbulent  boundary  layer  interactions  with  significant  suc- 
cess. Examples  include  Adams,  1998,  Urbin  et  al.,  1999,  Rizzetta  et  ah, 
2000  and  Rizzetta  and  Visbal,  2001. 

The  objective  of  this  paper  is  to  assess  the  capability  of  our  LES 
methodology  to  accurately  predict  the  flowfield  in  a supersonic  expansion- 
compression  corner  (Fig.  1).  This  configuration  is  reminiscent  of  aerody- 
namic configurations  wherein  a supersonic  boundary  layer  is  subjected 
to  an  initial  expansion  followed  by  a subsequent  compression.  Interest 
in  this  configuration  is  due  in  part  to  the  stabilizing  influence  of  the 
expansion  (Dussauge  and  Gaviglio,  1987;  Zheltovodov  et  al.,  1987;  Zhel- 
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tovodov  and  Schuelein,  1988;  Smith  and  Smits,  1997;  Stephen  et  al., 
1998;  Zheltovodov  et  al.,  1990b).  The  first  systematic  combined  ex- 
perimental and  numerical  study  of  an  expansion-compression  corner  by 
Zheltovodov  et  al.,  1992  and  Zheltovodov  et  al.,  1993  showed  that  several 
different  turbulence  models  (including  k — e,  q—u>  and  several  modifica- 
tions thereto)  did  not  accurately  predict  the  separation  and  attachment 
positions,  and  distributions  of  surface  skin  friction  and  heat  transfer.  We 
therefore  seek  to  ascertain  the  capability  of  LES  to  predict  this  flowfield. 


Moo 


Governing  Equations 

The  governing  equations  are  the  spatially  filtered,  Favre-averaged 
compressible  Navier-Stokes  equations.  The  spatial  filtering  removes  the 
small  scale  ( subgrid  scale)  fluctuations,  while  the  three  dimensional,  time 
dependent  large  scale  ( resolved  scale ) motion  is  retained.  For  an  arbi- 
trary function  T{xi,  /),  the  ordinary  and  Favre-filtered  variables  T{xt)  t) 
and  T{xi,  t ) are 


F{x, 


,t)=  [ G(xi 
Jd 


A).F(&,t)d&  and  T(xi,t) 


PT_ 

P 


(1) 


where  G is  the  filter  function,  and  A is  a measure  of  the  filter  width  and 
is  related  to  the  computational  mesh  size. 

The  filtered  governing  equations  using  Cartesian  tensor  notation  are 


dpuj 

dt 

dpe 

dt 


+ 


dp  dpui 
dt  dxi 

(2) 

dpuiUj  dp  i dTij 

dxj  dxi  dxj 

(3) 

d , _ dHj 

dxj(Pe+p)uj=  gxj 

(4) 

p - pRT 

(5) 
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where  xt  represents  the  Cartesian  coordinates  (i  = 1, 2, 3),  p is  the  mean 
density,  Hi  are  the  Cartesian  components  of  the  filtered  velocity  and  p 
is  the  mean  pressure.  The  total  stress  tensor  Tt]  = tvj  + o jj  where  the 
Subgrid  Scale  (SGS)  stress  n3  and  viscous  stress  al3  are 


Tij 


= -p(UiUj  - UiUj) 


P(T) 


2 duk 


duj 

dx-; 


3 Wj*  + ~ + — 


(6) 


where  p(T)  is  the  molecular  viscosity.  The  sum  of  the  heat  flux  plus 
work  done  by  the  stresses  is  Hj  = Qj  + q3  + 7 IjUi  where  the  SGS  and 
molecular  heat  fluxes  are 

Qj  = - cpp  (ujT  - ujf)  and  q3  = k(T)-^-  (7) 

where  k(T)  the  molecular  thermal  conductivity.  The  form  of  TLj  was 
proposed  by  Knight  et  al.,  1998  and  found  to  provide  an  accurate  model 
of  SGS  turbulent  diffusion  in  decaying  compressible  isotropic  turbulence 
(Martin  et  al.,  1999).  The  total  energy  pe  and  SGS  turbulence  kinetic 
energy  pk  per  unit  volume  are 

pe  = pcvf  + \pUiUi  + pk  and  pk  = \p  {vQui  - UiUi ) (8) 

Closure  of  the  system  of  equations  (2)  to  (5)  requires  specification  of 
a model  for  the  subgrid  scale  stress  and  heat  flux  Qj.  There  are  two 
basic  approaches  (Ghosal,  1999),  namely,  1)  the  explicit  specification  of 
an  SGS  model,  and  2)  the  Monotone  Integrated  Large  Eddy  Simulation 
(MILES)  method.  In  first  approach,  an  explicit  mathematical  model  for 
Tij  and  heat  flux  Q3  is  defined  ( e.g .,  the  Smagorinsky  model).  Examples 
are  presented  in  the  recent  reviews  of  Galperin  and  Orszag,  1993  and 
Lesieur  and  Metais,  1996.  In  the  second  approach,  the  SGS  model  is 
inherent  in  the  numerical  algorithm  (Boris  et  al.,  1992;  Oran  and  Boris, 
1993;  Grinstein,  1996;  Grinstein  and  Fureby,  1998;  Fureby  and  Grinstein, 
2000).  Fureby  and  Grinstein,  1999  showed  that  MILES  introduces  a 
tensor  eddy  diffusivity  into  the  equivalent  SGS  stress,  in  contrast  to  the 
isotropic  eddy  diffusivity  of  the  standard  explicit  Smagorinsky-type  SGS 
models. 

The  no-slip  condition  is  applied  at  solid  (impermeable)  boundaries. 
The  downstream  boundary  condition  for  supersonic  flows  is  typically  a 
zero  gradient  condition  on  the  conservative  flow  variables  (p,  put , pe). 
Periodic  boundary  conditions  are  usually  employed  for  the  spanwise 
boundaries  with  the  requirement  that  the  spanwise  domain  is  large  com- 
pared to  the  energy  containing  eddies  of  the  flow  and  the  flowfield  is  sta- 
tistically homogeneous  in  the  spanwise  direction.  The  farfield  boundary 
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condition  for  supersonic  flows  is  typically  a Riemann  condition  allowing 
waves  to  leave  the  computational  domain  without  reflection.  The  inflow 
boundary  condition  for  boundary  layers  is  a time-dependent  boundary 
layer  profile  obtained  by  the  rescaling  method  originally  developed  by 
Lund  et  al.,  1998  for  incompressible  boundary  layers  and  extended  to 
compressible  boundary  layers  by  Urbin  and  Knight,  1999.  The  initial 
condition  is  typically  obtained  by  linear  interpolation  from  a previous 
simulation  at  comparable  Mach  and  Reynolds  numbers. 

Numerical  Algorithm 

The  governing  equations  (2)  to  (5)  are  solved  using  a unstructured 
grid  of  tetrahedra.  The  finite  volume  algorithm  is  second  order  accurate 
in  space  and  time.  The  inviscid  fluxes  are  computed  using  Godunov’s 
method  with  the  left  and  right  states  at  each  face  reconstructed  using 
a second  order  Least  Squares  method  (Okong’o  and  Knight,  1998a). 
The  stencil  of  cells  employed  for  reconstruction  is  isotropic  except  in 
the  vicinity  of  shock  waves  where  an  ENO-like  anisotropic  stencil  is 
employed  (Chernyavsky  et  al.,  2001).  The  MILES  methodology  is  em- 
ployed (i.e.,  Tij  = 0,  Qj  — 0).  The  molecular  viscous  stresses  and  heat 
flux  are  obtained  using  a discrete  version  of  Gauss’  Theorem  (Okong’o 
and  Knight,  1998a).  The  temporal  integration  is  performed  by  using 
a second-order  accurate  Runge-Kutta  method.  The  code  is  parallelized 
using  domain  decomposition  with  the  Message  Passing  Interface  (MPI). 
The  flow  variables  are  non-dimensionalized  using  the  incoming  bound- 
ary layer  thickness  S,  and  incoming  freestream  velocity  Uoo,  density  px, 
static  temperature  T and  molecular  viscosity  fi^. 

The  code  has  been  validated  for  a variety  of  turbulent  flows  by  com- 
parison with  experiment  and  Direct  Numerical  Simulation  (DNS).  Ex- 
amples include  decay  of  isotropic  turbulence  (Knight  et  al.,  1997;  Knight 
et  al.,  1998),  incompressible  channel  flow  (Okong’o  and  Knight,  1998b; 
Okong’o  et  al.,  2000),  supersonic  turbulent  boundary  layer  (Urbin  et  al., 
1999;  Urbin  and  Knight,  1999;  Yan  et  al.,  2000;  Urbin  and  Knight, 
2001),  and  supersonic  compression  corner  (Urbin  et  al.,  1999;  Urbin 
et  al.,  2000;  Yan  et  al.,  2000;  Chernyavsky  et  al.,  2001).  The  supersonic 
boundary  layer  results  are  summarized  in  the  next  section. 

Flat  Plate  Boundary  Layer 

Urbin  and  Knight,  2001  performed  an  LES  of  an  adiabatic  Mach  3 
boundary  layer.  A detailed  grid  refinement  study  was  performed  to  as- 
certain the  required  grid  resolution  in  the  viscous  sublayer,  logarithmic 
and  outer  regions  of  the  boundary  layer.  The  computed  mean  veloc- 
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ity  profile,  expressed  in  Van  Driest  transformed  notation,  is  shown  in 
Fig.  2.  The  profile  shows  excellent  agreement  with  the  logarithmic  re- 
gion of  the  Law  of  the  Wall.  The  computed  adiabatic  wall  temperature 
is  within  3%  of  the  empirical  formula  Taw  = ^ (l  + i(7  - 1 )PrtmMl) 
where  Prtm  = 0.89  is  the  mean  turbulent  Prandtl  number.  The  com- 
puted friction  velocity  uT  is  within  5%  of  the  correlation  obtained  from 
the  combined  Law  of  the  Wall  and  Wake.  The  computed  normalized 
Reynolds  shear  stress  < p ><  u"v"  > /tw,  shown  in  Fig.  3,  shows 
excellent  agreement  with  the  experimental  data. 


y/8 


Figure  2.  Mean  Van  Driest  velocity 


Figure  3.  Reynolds  shear  stress 


Expansion-Compression  Corner 

Details  of  Computation 

The  flowfield  configuration  is  shown  in  Fig.  1.  An  incoming  Mach  3 
adiabatic  equilibrium  turbulent  boundary  layer  of  height  5 expands  over 
a 25°  corner  followed  by  a 25°  compression.  The  distance  along  the  ex- 
pansion surface  is  7.15  ( i.e .,  the  vertical  distance  between  the  two  hori- 
zontal surfaces  is  35,  and  the  horizontal  distance  between  the  expansion 
and  compression  corners  is  6.435). 

The  Cartesian  coordinates  x,y  and  z are  aligned  in  the  incoming 
streamwise,  transverse  and  spanwise  directions  with  the  origin  at  the 
inflow  boundary.  The  computational  domain  is  Lx  = 24.05,  Ly  = 3.45, 
and  Lz  = 1.9255.  The  expansion  corner  is  located  at  45  from  the  in- 
flow boundary.  The  grid  consists  of  253  x 35  x 57  nodes  in  the  x,  y and 
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z directions,  respectively,  forming  479,808  hexahedra  which  are  subdi- 
vided into  five  tetrahedra  each.  Thus,  the  total  number  of  tetrahedra  is 
2,399,040.  The  grid  is  stretched  in  the  y direction  with  spacing  0.008<5 
at  the  wall  and  a geometric  stretching  factor  of  1.154.  The  grid  is  con- 
centrated in  the  streamwise  direction  in  the  neighborhood  of  the  expan- 
sion and  compression  corners.  The  details  are  shown  in  Table  1 where 
A y+  = AyuT/vw  where  vw  is  the  computed  kinematic  viscosity  at  the 
wall,  uT  = \J tw  / pw  is  the  friction  velocity,  tw  is  the  computed  wall  shear 
stress  and  pw  is  the  computed  density  at  the  wall.  The  grid  is  consistent 
with  the  resolution  requirements  for  the  LES  code  established  by  Urbin 
and  Knight,  2001. 

Table  1.  Details  of  Grid 


Name 

Ax+ 

A y+ 

at  the  wall 

A z+ 

Ax/S 

Ay/5 
at  y = S 

Az/S 

Tetras 

Computed 

20.9 

1.67 

D 

0.1 

0.14 

0.034 

2,399,040 

The  inflow  boundary  condition  is  obtained  from  a separate  flat  plate 
boundary  layer  computation.  All  the  quantities  are  averaged  in  time 
and  in  the  spanwise  direction  and  denoted  by  < / >.  The  time  av- 
eraging period  is  set  to  three  times  the  flow-through  time,  where  one 
flow-through  time  is  defined  as  the  time  for  the  freestream  flow  to  tra- 
verse the  computational  domain.  The  averaging  is  performed  once  the 
initial  transient  has  decayed  (i.e.,  after  four  flow-through  times).  The 
details  are  presented  in  Urbin  et  al.,  1999. 

Experiments 

Experimental  data  has  been  obtained  by  Zheltovodov  et  al.,  1987, 
Zheltovodov  and  Schuelein,  1988,  and  Zheltovodov  et  al.,  1990a  and 
presented  in  part  in  tabular  form  in  Zheltovodov  et  al.,  1990b  for  the 
expansion-compression  corner  at  Mach  3 and  several  Reynolds  numbers 
Reg  based  on  the  incoming  boundary  layer  thickness  6.  The  experimen- 
tal conditions  are  listed  in  Table  2,  where  FPBL  and  ECC  imply  flat 
plate  boundary  layer  and  expansion-compression  corner,  respectively. 
The  LES  was  performed  at  a lower  Reynolds  number  (Reg  = 2 x 104) 
than  the  experiment  (Reg  = 4.4  x 104  to  1.94  x 105)  for  reasons  of 
computational  cost.  Additional  LES  cases  will  be  performed  at  higher 
Reynolds  numbers. 
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Table  2.  Details  of  Experiments  and  Computation 


Cases 

Mach 

Res 

References 

ECC 

2.9 

4.07  x 104 

Zheltovodov  et  al. , 1990a 

ECC 

2.9 

6.76  x 104 

Zheltovodov  et  al.,  1990a 

ECC 

2.9 

8.0  x 104 

Zheltovodov  et  al.,  1990a 

ECC 

2.9 

1.94  x 105 

Zheltovodov  et  al.,  1987;  Zheltovodov  et  al.,  1990b 

ECC 

2.88 

2.0  x 104 

Present  computation 

FPBL 

2.88 

1.33  x 105 

Zheltovodov  et  al.,  1990b 

Results 


The  structure  of  the  flowfield  is  shown  in  Figs.  4 and  Fig.  5 which 
display  the  mean  static  pressure  and  streamlines  at  z = 6.  The  flow 
expands  around  the  first,  corner,  and  recompresses  at  the  second  cor- 
ner through  a shock  which  separates  the  boundary  layer  as  evident  in 
Fig.  5.  The  flowfield  structure  is  in  good  agreement  with  the  results  of 
Zheltovodov  et  ah,  1987;  Zheltovodov  and  Schuelein,  1988;  Zheltovodov 
et  ah,  1990a  and  Zheltovodov  et  ah,  1990b  which  are  shown  qualitatively 
in  Fig.  1. 


p/p. 


z=1.08 


9 0.94 

8 0.86 
7 0.78 

6 0.68 
5 0.60 

4 0.52 

3 0.44 

2 0.35 

Shock  wave  1 0.27 


2 


0 

-2 


to 

>< 


x/5 


Figure  4 • Mean  static  pressure  (s  is  separation,  a is  attachment) 
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Figure  5.  Mean  streamlines  (s  is  separation,  A is  attachment) 


The  mean  velocity  profiles  in  the  ^-direction  are  shown  in  Fig.  6 at 
x = 28  and  x = 65,  where  x is  measured  from  the  inflow  along  the 
direction  of  the  inflow  freestream  velocity  (Fig.  4).  The  abscissa  is  the 
component  of  velocity  locally  parallel  to  the  wall,  and  the  ordinate  is 
the  distance  measured  normal  to  the  wall.  The  first  profile  is  upstream 
of  the  expansion  corner  which  is  located  at  x = 45,  and  the  second  is 
downstream  of  the  expansion  fan  and  upstream  of  the  separation  point. 
The  computed  mean  velocity  profile  at  the  first  location  is  slightly  fuller 
than  the  experiment.  This  is  consistent  with  the  experimentally  observed 
dependence  of  the  exponent  n in  the  power- law  U/Uoo  = (y/5)1/71  on  the 
Reynolds  number.  The  second  profile  shows  a significant  acceleration  of 
the  flow  in  the  outer  portion  of  the  boundary  layer  due  to  the  expansion. 


Figure  6.  Mean  velocity  Figure  7.  Separation  length 


<p»<u’u’>/p_U* 
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Zheltovodov  and  Schuelein,  1988  and  Zheltovodov  et  al.,  1993  devel- 
oped an  empirical  correlation  for  the  separation  length  (defined  as  the 
minimum  distance  between  the  mean  separation  and  attachment  points 
on  the  wall)  in  the  expansion-compression  corner  interaction.  The  scaled 
separation  length  Lsep/Lc  is  observed  experimentally  to  be  a function  of 
Re5  where  the  characteristic  length  (Lc)  is  defined  by 

Lc  = 6e(P2/ppl)31/M*  (9) 

where  5e  is  the  incoming  boundary  layer  thickness  (upstream  of  the 
expansion  corner),  p2  is  the  pressure  after  the  shock  in  inviscid  flow, 
ppl  is  the  plateau  pressure  from  the  empirical  formula ppi  = pe(^Me  + 1) 
where  pe  and  Me  are  the  static  pressure  and  freestream  Mach  number 
upstream  of  the  compression  corner  and  downstream  of  the  expansion 
fan.  In  the  computation,  the  location  is  taken  to  be  x = 65.  The 
values  of  Me  and  p2  have  been  computed  using  inviscid  theory.  Also, 
Rese  = 1.8  x 104  for  LES  ( Rese  = pet/e5e/jue,  where  pe,Ue  and  pe  are 
computed  using  inviscid  theory).  The  experimental  data  correlation  of 
Zheltovodov  and  Schuelein,  1988  and  the  computed  result1  for  the  scaled 
separation  length  is  shown  in  Fig.  7.  The  computed  value  is  consistent 
with  a linear  extrapolation  of  the  expermental  data. 


Figure  8.  Reynolds  streamwise  stress  Figure  9.  Surface  pressure 

The  surface  pressure  profile  in  Fig.  9 displays  a pressure  plateau  on 
the  compression  face  generated  by  the  separation  bubble.  The  exper- 


^he  uncertainty  in  the  computed  value  of  Lsep/Lc  is  associated  with  the  uncertainty  in 
determining  8e.  We  have  used  the  streamwise  Reynolds  stress  (<  p ><  u'u'  >)  to  determine 
6e  (Fig.  8),  where  u'  is  the  fluctuating  velocity  parallel  to  the  wall. 
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iments  exhibit  a trend  of  increase  in  the  size  of  the  pressure  plateau 
region  with  decreasing  Reynolds  number.  The  experimental  data  at  the 
lowest  Reynolds  number  (Res  =4.1  x 104)  shows  close  agreement  with 
the  computed  results  for  Reg  = 2 x 104  for  the  location,  extent  and  mag- 
nitude of  the  pressure  plateau.  Moreover,  the  shape  of  the  experimental 
pressure  plateau  shows  little  variation  for  Reg  < 6.8  x 104,  thus  suggest- 
ing that  the  computed  pressure  plateau  region  (for  Reg  = 2 x 104)  is 
accurate.  The  computed  recovery  of  the  surface  pressure  is  more  rapid 
than  in  the  experiment,  however. 


The  computed  and  experimental 
mean  skin  friction  coefficient  Cf  = 
Tw/^PooU^  are  shown  in  Fig.  10. 
The  computed  separation  and  at- 
tachment points  are  evident.  The 
skin  friction  rises  rapidly  down- 
stream of  attachment.  The  com- 
puted results  at  Reg  = 2 x 104  are 
in  close  agreement  with  the  experi- 
mental data  at  Reg  = 8.0  x 104  and 
1.94  x 105  in  the  region  downstream 
of  reattachment. 


Figure  10.  Skin  friction  coefficient 


Summary 

An  unstructured  grid  Large  Eddy  Simulation  methodology  has  been 
validated  for  complex  compressible  turbulent  flows.  The  methodology 
is  based  on  the  MILES  concept  wherein  the  inherent  dissipation  of  the 
monotone  inviscid  flux  algorithm  provides  the  energy  transfer  from  the 
resolved  to  the  subgrid  scales.  The  methodology  has  been  validated  by 
comparison  with  experiment  for  a variety  of  supersonic  turbulent  flows 
including  a turbulent  boundary  layer  and  an  expansion-compression  cor- 
ner at  Mach  3. 
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